In [1]:
import plotly
import plotly.graph_objs as go
import plotly.express as px
from plotly.subplots import make_subplots

import torch
import numpy as np
import pandas as pd

import dill as pickle

import random
import math

from tedeous.solver import grid_format_prepare
import statistics
import scipy.io
import json
In [2]:
def get_rms(records):
    """
        Root-mean-square (rms)
    """
    return math.sqrt(sum([x ** 2 for x in records]) / len(records))
In [3]:
path = ""
mat = scipy.io.loadmat(f'{path}burgers.mat')

data = mat['u']
data = np.transpose(data)
t = np.ravel(mat['t'])
x = np.ravel(mat['x'])
In [4]:
t_torch = torch.from_numpy(t).float()
x_torch = torch.from_numpy(x).float()

param = [x_torch, t_torch]
In [5]:
prepared_grid_main = grid_format_prepare(param, "mat")
In [6]:
data.shape
Out[6]:
(256, 256)
In [7]:
fig = go.Figure(data=[go.Surface(x=x, y=t, z=data, name='Initial field', legendgroup='i',
                                 opacity=0.5)])
    
fig.update_layout(scene_aspectmode='auto')
fig.update_layout(scene=dict(
                      xaxis_title='x1 = x',
                      yaxis_title='x2 = t',
                      zaxis_title='u',
                      aspectratio={"x": 1, "y": 1, "z": 1}
                  ))
fig.show()
In [8]:
fixed_t_index = 0
fixed_t = t[fixed_t_index]
fixed_data_slice = data[fixed_t_index, :]

fig = go.Figure(data=[go.Scatter(x=x, y=fixed_data_slice)])

fig.update_layout(
    title=f"График данных при t = {fixed_t}",
    xaxis_title="X",
    yaxis_title="Значение данных",
)

fig.show()
In [9]:
set_solutions = torch.load(f'{path}file_u_main_[35, 256, 256]_mat_[0].pt')
In [10]:
u = data
In [11]:
for k in range(len(set_solutions)):
    error_rmse = np.sqrt(np.mean((u.reshape(-1) - set_solutions[k].reshape(-1)) ** 2))
    print(f'{k + 1}. RMSE = {error_rmse}')

1. RMSE = 82.25906457258895
2. RMSE = 82.25906457258895
3. RMSE = 82.25906457258895
4. RMSE = 82.25906457258895
5. RMSE = 82.25906457258895
6. RMSE = 82.25906457258895
7. RMSE = 82.21910825931126
8. RMSE = 82.21910825931126
9. RMSE = 82.21910825931126
10. RMSE = 82.21910825931126
11. RMSE = 82.21910825931126
12. RMSE = 82.21910825931126
13. RMSE = 82.42049484571486
14. RMSE = 82.21910825931126
15. RMSE = 82.21910825931126
16. RMSE = 82.42049484571486
17. RMSE = 82.21910825931126
18. RMSE = 82.21910825931126
19. RMSE = 82.42049484571486
20. RMSE = 82.21910825931126
21. RMSE = 82.42049484571486
22. RMSE = 82.21910825931126
23. RMSE = 82.21910825931126
24. RMSE = 82.21910825931126
25. RMSE = 82.21910825931126
26. RMSE = 82.42049484571486
27. RMSE = 82.42049484571486
28. RMSE = 82.21910825931126
29. RMSE = 82.42049484571486
30. RMSE = 82.42049484571486
31. RMSE = 82.21910825931126
32. RMSE = 82.21910825931126
33. RMSE = 82.21910825931126
34. RMSE = 82.21910825931126
35. RMSE = 82.21910825931126
In [12]:
mean_arr = np.zeros((set_solutions.shape[1], set_solutions.shape[2]))
var_arr = np.zeros((set_solutions.shape[1], set_solutions.shape[2]))
s_arr = np.zeros((set_solutions.shape[1], set_solutions.shape[2]))  # sample standard deviation of data
In [13]:
mean_arr = np.zeros((set_solutions.shape[1], set_solutions.shape[2]))
var_arr = np.zeros((set_solutions.shape[1], set_solutions.shape[2]))
s_arr = np.zeros((set_solutions.shape[1], set_solutions.shape[2]))  # sample standard deviation of data
In [14]:
for i in range(set_solutions.shape[1]):
    for j in range(set_solutions.shape[2]):
        mean_arr[i, j] = np.mean(set_solutions[:, i, j])
        var_arr[i, j] = np.var(set_solutions[:, i, j])
        s_arr[i, j] = np.std(set_solutions[:, i, j], ddof=1) #statistics.stdev()
In [15]:
mean_tens = torch.from_numpy(mean_arr)
var_tens = torch.from_numpy(var_arr)
s_arr = torch.from_numpy(s_arr)
In [16]:
# Confidence region for the mean
upper_bound = mean_tens + 1.96 * s_arr / math.sqrt(len(set_solutions))
lower_bound = mean_tens - 1.96 * s_arr / math.sqrt(len(set_solutions))
In [17]:
confidence_region = torch.cat((upper_bound, torch.flip(lower_bound, dims=(0,))), 0)
confidence_grid = torch.cat((prepared_grid_main.reshape(-1),  torch.flip(prepared_grid_main.reshape(-1), dims=(0,))), 0)
In [18]:
mean_tens = mean_tens.reshape(-1)
upper_bound = upper_bound.reshape(-1)
lower_bound = lower_bound.reshape(-1)
In [19]:
confidence_grid.shape
Out[19]:
torch.Size([262144])
In [20]:
fig = go.Figure(data=[
    go.Mesh3d(x=prepared_grid_main[0].reshape(-1), y=prepared_grid_main[1].reshape(-1), z=mean_tens, name='Solution field',
              legendgroup='s', showlegend=True, color='lightpink',
              opacity=1),
    go.Mesh3d(x=prepared_grid_main[0].reshape(-1), y=prepared_grid_main[1].reshape(-1), z=upper_bound, name='Confidence region',
              legendgroup='c', showlegend=True, color='blue',
              opacity=0.20),
    go.Mesh3d(x=prepared_grid_main[0].reshape(-1), y=prepared_grid_main[1].reshape(-1), z=lower_bound, name='Confidence region',
              legendgroup='c', color='blue', opacity=0.20)
])

fig.add_trace(go.Mesh3d(x=prepared_grid_main[0].reshape(-1), y=prepared_grid_main[1].reshape(-1), z=torch.from_numpy(u).reshape(-1),
              name='Initial field',
              legendgroup='i', showlegend=True, color='rgb(139,224,164)',
              opacity=0.5))

fig.update_layout(scene_aspectmode='auto')
fig.update_layout(showlegend=True,
                  scene=dict(
                      xaxis_title='x1',
                      yaxis_title='x2',
                      zaxis_title='u',
                      aspectratio={"x": 1, "y": 1, "z": 1}
                  ),
                  height=800, width=800
                  )

fig.show()
In [25]:
fixed_t_index = 100
fixed_t = t[fixed_t_index]
fixed_mean_slice = mean_tens.reshape(data.shape)[fixed_t_index, :]
fixed_upper_slice = upper_bound.reshape(data.shape)[fixed_t_index, :]
fixed_lower_slice = lower_bound.reshape(data.shape)[fixed_t_index, :]

fig = go.Figure(data=[go.Scatter(x=x, y=fixed_mean_slice), go.Scatter(x=x, y=data[fixed_t_index, :]),
                     go.Scatter(x=x, y=fixed_upper_slice), go.Scatter(x=x, y=fixed_lower_slice)])

fig.update_layout(
    title=f"График данных при t = {fixed_t}",
    xaxis_title="X",
    yaxis_title="Значение данных",
)

fig.show()
In [27]:
fig = make_subplots(rows=2, cols=1) #, subplot_titles=("Среднее решение", "Исходные данные"))


fig.add_trace(
    go.Contour(x=prepared_grid_main[0].reshape(-1),
                           y=prepared_grid_main[1].reshape(-1),
                           z=mean_tens, zmin=0, zmax=1000,
                           contours_coloring='heatmap',
                            showscale=False), row=1, col=1
)


# 'Визуализация "среднего" решения'
fig.add_trace(
    go.Contour(x=prepared_grid_main[0].reshape(-1),
                           y=prepared_grid_main[1].reshape(-1),
                           z=data.reshape(-1),
                           colorbar_len=1.04, # 1.06
                           contours_coloring='heatmap'), row=2, col=1
)

# 'Визуализация исходных данных'
fig.update_layout(paper_bgcolor='rgba(0,0,0,0)',
                  plot_bgcolor='rgba(0,0,0,0)',
                  margin=dict(l=0, r=0, t=0, b=0),
                  legend=dict(yanchor="top", y=0.98, xanchor="right", x=0.99, bgcolor = 'white'))

fig.update_xaxes(title_text="x", row=1, col=1)
fig.update_xaxes(title_text="x", row=2, col=1)

# Update yaxis properties
fig.update_yaxes(title_text="t", row=1, col=1)
fig.update_yaxes(title_text="t", row=2, col=1)

fig.update_layout(
    plot_bgcolor='white',
#     title_text='Объединенные графики',
    showlegend=True,
)
fig.update_xaxes(
    mirror=True,
    ticks='outside',
    showline=True,
    linecolor='black',
    gridcolor='lightgrey'
)
fig.update_yaxes(
    mirror=True,
    ticks='outside',
    showline=True,
    linecolor='black',
    gridcolor='lightgrey'
)

    
fig.show()
In [28]:
# plotly.io.write_image(fig, 'output_file_final.pdf', format='pdf')